﻿#include <xmmintrin.h> //SSE
#include <emmintrin.h> //SSE2
#include <pmmintrin.h> //SSE3
#include <tmmintrin.h> //SSSE3
#include <smmintrin.h> //SSE4.1
#include <nmmintrin.h> //SSSE4.2
#include <immintrin.h> //AVX¡¢AVX2
#include <pthread.h>
#include <stdlib.h>
#include <iostream>
#include <sys/time.h>
#include <semaphore.h>
#include <unistd.h>
#include <omp.h>
using namespace std;

#define NUM_THREADS 7

void OpenMP_Guided_avx(float** a, int n)
{
#pragma omp parallel for simd num_threads(NUM_THREADS) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_Guided(float** a, int n)
{
#pragma omp parallel for num_threads(NUM_THREADS) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_Dynamic_avx(float** a, int n)
{
#pragma omp parallel for simd num_threads(NUM_THREADS) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_Dynamic(float** a, int n)
{
#pragma omp parallel for num_threads(NUM_THREADS) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_Static_avx(float** a, int n)
{
#pragma omp parallel for simd num_threads(NUM_THREADS) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_Static(float** a, int n)
{
#pragma omp parallel for num_threads(NUM_THREADS) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_avx(float** a, int n)
{
#pragma omp simd
    for (int k = 0; k < n; k++)
    {

        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
        for (int i = k + 1; i < n; i++)
            a[i][k] = 0;
    }
}
void Serial(float** a, int n)
{
    int i, j, k;
    for (k = 0; k < n; k++)
    {
        for (j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (i = k + 1; i < n; i++)
        {
            for (j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void avx(float** a, int n)
{
    int i, j, k;
    for (k = 0; k < n; k++)
    {
        for (j = k + 1; j % 8 != 0 && j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        for (; j <= n - 8; j += 8)
        {
            __m256 temp;
            __m256 t = _mm256_set1_ps(a[k][k]);
            temp = _mm256_loadu_ps(a[k] + j);
            temp = _mm256_div_ps(temp, t);
            _mm256_storeu_ps(a[k] + j, temp);
        }
        for (; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (i = k + 1; i < n; i++)
        {
            for (j = k + 1; j % 8 != 0 && j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            for (; j <= n - 8; j += 8)
            {
                __m256 t1;
                __m256 t2;
                t1 = _mm256_loadu_ps(a[i] + j);
                t2 = _mm256_loadu_ps(a[k] + j);
                __m256 t3 = _mm256_set1_ps(a[i][k]);
                t3 = _mm256_mul_ps(t3, t2);
                t1 = _mm256_sub_ps(t1, t3);
                _mm256_storeu_ps(a[i] + j, t1);
            }
            for (; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_Guided_sse(float** a, int n)
{
#pragma omp parallel for simd num_threads(NUM_THREADS) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_Dynamic_sse(float** a, int n)
{
#pragma omp parallel for simd num_threads(NUM_THREADS) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_Static_sse(float** a, int n)
{
#pragma omp parallel for simd num_threads(NUM_THREADS) shared(a,n)
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_sse(float** a, int n)
{
#pragma omp simd
    for (int k = 0; k < n; k++)
    {

        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
        for (int i = k + 1; i < n; i++)
            a[i][k] = 0;
    }
}
void sse(float** a, int n)
{
    int i, j, k;
    for (k = 0; k < n; k++)
    {
        for (j = k + 1; j % 4 != 0 && j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        for (; j <= n - 4; j += 4)
        {
            __m128 temp;
            __m128 t = _mm_set1_ps(a[k][k]);
            temp = _mm_loadu_ps(a[k] + j);
            temp = _mm_div_ps(temp, t);
            _mm_storeu_ps(a[k] + j, temp);
        }
        for (; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (i = k + 1; i < n; i++)
        {
            for (j = k + 1; j % 4 != 0 && j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            for (; j <= n - 4; j += 4)
            {
                __m128 t1;
                __m128 t2;
                t1 = _mm_loadu_ps(a[i] + j);
                t2 = _mm_loadu_ps(a[k] + j);
                __m128 t3 = _mm_set1_ps(a[i][k]);
                t3 = _mm_mul_ps(t3, t2);
                t1 = _mm_sub_ps(t1, t3);
                _mm_storeu_ps(a[i] + j, t1);
            }
            for (; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void m_rest(float** a, int n)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < i; j++)
            a[i][j] = 0;
        a[i][i] = 1.0;
        for (int j = i + 1; j < i; j++)
        {
            float x = rand();
            while (x == 0) x = rand();
            a[i][j] = x;
        }
    }
    for (int k = 0; k < n; k++)
        for (int i = k + 1; i < n; i++)
            for (int j = 0; j < n; j++)
                a[i][j] += a[k][j];
}
int main()
{
    srand((unsigned)time(NULL));
    int n = 10;
    while (n < 1600)
    {
        float** a = (float**)aligned_alloc(32, n * sizeof(float*));
        for (int i = 0; i < n; i++)
            a[i] = (float*)aligned_alloc(32, n * sizeof(float));
        m_rest(a, n);

        cout << n << " ";
        struct timeval tstart, tend;
        double timeUsed = 0.0;
        float** t = (float**)aligned_alloc(32, n * sizeof(float*));
        for (int i = 0; i < n; i++)
            t[i] = (float*)aligned_alloc(32, n * sizeof(float));

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        Serial(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        double time1 = timeUsed;
        cout << timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Static_avx(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1 / timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Dynamic_avx(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1 / timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Guided_avx(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1 / timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Static_sse(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1 / timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Dynamic_sse(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1 / timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Guided_sse(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1 / timeUsed << " ";

        cout << endl;
        for (int i = 0; i < n; i++)
        {
            delete[] t[i];
            delete[] a[i];
        }
        delete[] t;
        delete[] a;
        if (n >= 800) n += 200;
        if (n >= 400 && n < 800) n += 100;
        if (n >= 200 && n < 400) n += 50;
        if (n < 200) n += 10;
    }
    return 0;
}